%{
load variable from ISIMIP2b dataset, take LSWT as an example
%}

var = { 'surftemp_mean' };

for flag_rcp = [1:5]      %    #1 = historical       #2 = picontrol        #3 = rcp26      #4 = rcp60      #5 = rcp85

    disp( rcp{flag_rcp} )

    for flag_period = [1:5]   %    %    #1 = DJF       #2 = MAM       #3 = JJA      #4 = SON      #5 = annual

        for flag_model = [1:3]      %   #1 = albm     #2 = lake      #3 = simstrat-uog

            for flag_gcm = [1:4]      %   #1 = gfdl-esm2m      #2 = hadgem2-es      #3 = ipsl-cm5a-lr      #4 = miroc5

                %%% file name for rcp data
                loc_rcp = [ '/home/work_file/', ...
                    model{flag_model}, '/' ];
                file_rcp = [ period{ flag_period }, '_', model{ flag_model }, '_' gcm{flag_gcm}, '_', ...
                    rcp{ flag_rcp }, '_LSWT.mat' ];

                %%% load rcp data
                LSWT_rcp = load( [loc_rcp, file_rcp], var{:}) ;

                %%% save to ensemble members
                for flag_var = 1:length(var)

                    if flag_rcp == 1        % historical   

                        % joint historical and rcp60
                        file_rcp_60 = [ period{ flag_period }, '_', model{ flag_model }, '_' gcm{flag_gcm}, '_', ...
                                    rcp{ 4 }, '_LSWT.mat' ];
                        LSWT_rcp_60 = load( [loc_rcp, file_rcp_60], var{:}) ;
                        LSWT_rcp.( var{ flag_var } ) = cat( 3,LSWT_rcp.( var{ flag_var } ), LSWT_rcp_60.( var{ flag_var } ) );

                        LSWT_rcpMean_ensemble.(rcp{flag_rcp}).(period{flag_period}){ flag_model }{ flag_var }( :, :, flag_gcm ) = ... 
                                mean( LSWT_rcp.( var{ flag_var } )(:, :, time_present-1861+1 ), 3, "omitnan");
                    
                    elseif flag_rcp == 2        % picontrol

                        LSWT_rcpMean_ensemble.(rcp{flag_rcp}).(period{flag_period}){ flag_model }{ flag_var }( :, :, flag_gcm ) = ... 
                                mean( LSWT_rcp.( var{ flag_var } )(:, :, : ), 3, "omitnan");    % picontrol use the total time serise
                    
                    else        % RCPs

                        LSWT_rcpMean_ensemble.(rcp{flag_rcp}).(period{flag_period}){ flag_model }{ flag_var }( :, :, flag_gcm ) = ... 
                                mean( LSWT_rcp.( var{ flag_var } )(:, :, time_future-2006+1 ), 3, "omitnan");
                    
                    end
                    
                end

            end

        end

    end

end

clear loc_rcp file_rcp file_rcp_60 LSWT_rcp LSWT_rcp_60










